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Abstract 

New numerical methods have been applied in relativity to obtain a 
numerical evolution of Einstein equations much more robust and stable. 
Starting from 3+1 formalism and with the evolution equations written 
as a FOFCH (first-order flux conservative hyperbolic) system, advanced 
numerical methods from CFD (computational fluid dynamics) have been 
successfully applied. A flux limiter mechanism has been implemented in 
order to deal with steep gradients like the ones usually associated with 

■ black hole spacetimes. As a test bed, the method has been applied to 
' 3D metrics describing propagation of nonlinear gauge waves. Results are 

CN) , compared with the ones obtained with standard methods, showing a great 

■ increase in both robustness and stability of the numerical algorithm. 
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1 Introduction 



From the very beginning, 3D numerical relativity has not been an easy domain. 
Difficulties arise either from the computational side (the large amount of vari- 
. ables to evolve, the large number of operations to perform, the stability of the 

evolution code) or from the physical side, like the complexity of the Einstein 
equations themselves, boundary conditions, singularity avoidant gauge choices, 
and so on. Sometimes there is a connection between both sides. For instance, 
the use of singularity avoidant slicings generates large gradients in the vicinity 
of black holes. Numerical instabilities can be produced by these steep gradients. 
The reason for this is that the standard evolution algorithms are unable to deal 
with sharp profiles. The instability shows up in the form of spurious oscillations 
which usually grow and break down the code. 

Numerical advanced methods from CFD (Computational Fluid Dynamics) 
can be used to avoid this. Stable codes are obtained which evolve in a more 
robust way, without too much dissipation, so that the shape of the profiles of the 
evolved quantities is not lost. These advanced methods are then specially suited 
for the problem of shock propagation, but they apply only to strongly hyperbolic 
systems, where one is able get a full set of eigcnficlds which generates all the 
physical quantities to be evolved. In the ID case, these methods usually fulfill 
the TVD (Total Variation Diminishing) condition when applied to transport 
equations. This ensures that no new local extreme appear in the profiles of the 
eigenfields, so that spurious oscillations are ruled out ab initio (monotonicity 
preserving condition) . Unluckily, there is no general method with this property 
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in the 3D case, mainly because the eigenfield basis depends on the direction of 
propagation. We will show how this can be achieved at least in some cases. 

The specific methods we will use are known as flux limiter algorithms. We 
will consider plane waves in 3D as a first generalization of the ID case, be- 
cause the propagation direction is constant. This specific direction leads then 
to an specific eigenfield basis, so that the ID numerical method can be easily 
generalized to the 3D case. 

The algorithm will be checked with a "Minkowski waves" metric. It can be 
obtained by a coordinate transformation from Minkowski space-time. All the 
metric components are transported while preserving their initial profiles. The 
line element has the following form: 

ds 2 = -H(x - t) dt 2 + H(x - t) dx 2 + dy 2 + dz 2 (I) 

where H(x — t) is any positive function. We can choose a periodic profile with 
sharp peaks so both the space and the time derivatives of H(x — t) will have 
discontinuous step-like profiles. If we can solve well this case (the most extreme), 
we can hope that the algorithm will work as well in more realistic cases where 
discontinuities do not appear. 

Minkowski waves are a nice test bed because the instabilities can arise only 
from the gauge (there are pure gauge after all!). This is a first step to deal 
with evolution instabilities in the Einstein equations by the use of flux limiter 
methods. This will allow us to keep all our gauge freedom available to deal 
with more physical problems, like going to a co-rotating frame or adapting to 
some special geometry. Advanced numerical methods take care of numerical 
problems so that 'physical' gauge choices can be used to take care of physics 
requirements. 



2 The system of equations 

We will use the well known 3+1 description of spacetinie (j], 0, || which starts 
by decomposing the line element as follows: 

ds 2 = -a 2 dt 2 + 7y {dx 1 + I3 l dt) (dx> + (3 j dt) i,j = 1,2,3 (2) 

where 7y- is the metric induced on the three-dimensional slices and (3 t is the 
shift. For simplicity the case (3 l = (normal coordinates) will be considered. 
The intrinsic curvature of the slices is then given by the three-dimensional Ricci 
tensor ^.Ry, whereas their extrinsic curvature K%j is given by: 

dtHj = -2 a Kij (3) 

In what follows, all the geometrical operations (index raising, covariant deriva- 
tions, etc) will be performed in the framework of the intrinsic three-dimensional 
geometry of every constant time slice. With the help of the quantities defined 
in (Q,||) , the ten fourdimensional field equations can be expressed as a set of six 
evolution equations: 

dtKy = -ViOj + a [WRij - 2K 2 + tr K K tj - 8^ (T„- - | 7ii )] (4) 
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plus four constraint equations 



&R-tr(K 2 ) + (trK) 2 = IGira 2 T 00 (5) 

V fc K k i-di{trK) = MaT° i (6) 

The evolution system (||) has been used by numerical relativists since the very 
beginning of the field (see for instance the seminal work of Eardley and Smarr 
(U), both in spherically symmetric (ID) and axially symmetric (2D) spacetimes. 
By the turn of the century, the second order system (^) has been rewritten as a 
first-order flux conservative hyperbolic (FOFCH) system || ^, in order to deal 
with the generic (3D) case, where no symmetries are present. But the second 
order system (^) is still being used in 3D numerical calculations @] , mainly when 
combined with the conformal decomposition of as introduced by Shibata 
and Nakamura B [Tc[| . In the system(|]||) there is a degree of freedom to be 
fixed because the evolution equation for the lapse function a is not given. In 
the study of Black Holes, the slicing is usually chosen in order to avoid the 
singularity jljj : 

dt In a = —a Q (7) 

where: 

Q = f(a)trK (8) 

Three basic steps are needed to obtain a FOFCH system from the ADM system. 
First, one must introduce some new auxiliary variables to express the second 
order derivatives in space as first order. These new quantities correspond to the 
space derivatives: 

A k = d k In a , D kij = 1/2 d k ^ij (9) 
The evolution equations for these variables are: 

d t A k + d k (aftrK) = (10) 
d t D kij + 8 k (a Kij) = (11) 

At the second step the system is expressed in a first order balance law form 

d t u + d k F k (u) = S(u), (12) 

where the array u displays the set of independent variables to evolve and both 
"fluxes" F k and "sources" S are vector valued functions. At the third step 
another additional independent variable is introduced to obtain a strongly hy- 
perbolic system ftlf : 

V, = D lr r - D r „ (13) 

and its evolution equation is obtained using the definition of Kij from (Q) and 
switching space and time derivatives in the momentum constraint The 
result is an independent evolution equation for Vi while the previous definition 
( [l3| ) in terms of space derivatives can be instead be considered as a first integral 
of the extended system. The extended array u will then contain the following 
37 functions u = (a, 7^, K v] , A t , D ki j, Vi). 
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3 The numerical algorithm 



Due to the structure of the equations, the evolution(represented by the operator 



E(At)) described by (12) can be decomposed into two separate processes; the 
first one is a transport process and the second one is the contribution of the 
sources. 

The sources step (represented by the operator S(At)) does not involve space 
derivatives of the fields, so that it consists in a system of coupled non-linear 
ODE (Ordinary Differential Equations): 

d t u = S(u) (14) 

The transport step (represented by the operator T{At)) contains the principal 
part and it is given by a set of quasi- linear transport equations: 

d t u + d k F k {u) = (15) 

The numerical implementation of these separated processes is quite easy. Second 
order accuracy in At can be obtained by using the well known Strang splitting. 

E(At) = S(At/2) T(At) S(At/2) (16) 

According to (|[0) the lapse and the metric have no flux terms. It means that 
a reduced set of 30 quantities u = (Kij , Ai, D k ij , Vj) are transported in the 
second step over an inhomogeneous static background composed by (a, 7^). 
The equations for the transport step (fla) are given by: 



dtKij + d k (a\ k ij ) = (17) 

d t A k + d k (af{a)trK) = (18) 

d t D kij + d k {a K^) = (19) 

d t V k = (20) 

where: 

X\ 3 = D% - jV k ll3 + 1/2 5 k (Aj + 2Vj - D jr r ) + 1/2 5 k {A l + 2V> - D ir r ) (21) 

and m is an arbitrary parameter. 

To evolve the transport step, we will consider flux-conservative numeric al- 
gorithms p2] | , obtained by applying the balance to a single computational cell. 
In the ID case the cell goes from n to n+1 in time (t = n ■ At) and from j-1/2 
to j+1/2 in space (xj — j ■ Ax), so that we have: 

Interface fluxes can be calculated in many different ways, leading to different 
numerical methods. We will use here the well known MacCormack method. 
This flux-conservative standard algorithm works well for smooth profiles, as it 
can be appreciated in Figure 1. 

But this standard algorithm is not appropriate for step-like profiles because 
it produces spurious oscillations near the steep regions, as it can be appreciated 
in Figure 2. 
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Figure 1: Plot of K xx for the initial metric given by (|l|) with H{x — t) — 
1 + A cos[w (x — t))] with periodic boundaries. Continuous line is the initial 
condition. Dashed line is after 40 iterations 




Figure 2: Same as in Fig. 1 with the step- like initial data for K xx . Continuous 
line is the initial condition. Dashed line is after 10 iterations. Note the spurious 
oscillations around the corners 
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More advanced numerical methods must be used to eliminate (or at least 
to reduce) these oscillations. These advanced methods use information about 
the eigenfields and the propagation direction, so the flux characteristic matrix 
along the propagation direction must be diagonalized. 



4 Eigenfields 

We will use a convenient method to compute the eigenfields. Let us study the 
propagation of a step-like discontinuity in the transported variables u which will 
move along a specific direction n with a given velocity v. Information about 
the corresponding eigenfields can be extracted from the well known Rankine- 
Hugoniot shock conditions : 

v[u] = n k [F k (u)} (23) 

where [ ] represents the jump in the discontinuity. In our case 

v[Kij] = n r [a\ r l3 ] (24) 

v[A k ] = n k [af(a)trK] (25) 

v[D kij ] = n k [aK i:j ] (26) 

v[V k ] = (27) 

where we must note that both the background metric coefficients and the prop- 
agation direction are supposed to be continuous, so they are transparent to the 
] symbol. 

If we develop this expressions we arrive at the following conclusions, where 
S n = n r S r is the projection of the quantity S over n and S± — S k — S n n k are 
the transverse components: 

1) [V k ], [A±], [D±ij] and [A n — f trD n ] propagate along n with speed v = 0. 
There are 18 such eigenfields. For the line element given by (P n k is along the 
x axis and all these fields are actually zero. 

2) [\ n ij —trX n riinj] and [Kij — trKmrij] do generate eigenfields propagating 
along n with speed v = ±a (light cones). There are only 10 such eigenfields 
because all of them are traceless. For Minkowski waves, where there is only 
gauge, all these combinations are zero. This indicates that the correct way 
to get the traceless part of a given tensor Sij in this context is just to take 
Sij — trS rii Uj, so that the contribution of gauge modes will disappear. 

3) [A n ] and [trK] do generate eigenfields propagating along n with speed 
v = ±^/Ja (gauge cones). There are 2 such eigenfields corresponding to the 
gauge sector. For Minkowski waves, there are the only non-zero components. 
We are left with: 

v[trK] = a[A n ] (28) 
v[A k ] = af(a)[trK]n k (29) 

so that [A k ] is proportional to n k . Now we can get the gauge eigenfields: 

^fn k F k (trK)±F(A n ) (30) 

These eigenfields propagate along n according simple advection equations, a 
familiar situation in the ID case. Although this decomposition and diagonaliza- 
tion is trivial in ID, it is very useful in the multidimensional case for a generic 
direction n. 
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Figure 3: Same as Fig. 2 where the methods presented in this paper are applied. 
Continuous line is the initial condition. Dashed line is after 10 iterations 



5 Flux limiter methods 

The flux limiter methods jl^] we will use can be decomposed into some basic 
steps. First of all the interface fluxes have to be calculated with any standard 
second order accurate method (MacCormack in our case). Then, the propaga- 
tion direction n and the corresponding eigenfields can be properly identified at 
every cell interface. Two advection equations (one for every sense of propaga- 
tion) arc now available for the gauge eigenfluxes (|30|). 

Let us choose for instance the eigenflux which propagates to the right (an 
equivalent process will be valid for the other eigenflux propagating to the left). 

This interface eigenflux F™^j 2 can be understood as the grid point flux F" plus 

some increment Aj +1 / 2 = F™^/^ — F™. In general, the purpose of the limiter is 
to use of a mixture of the increments A J+1 / 2 and Aj_ 1 / 2 to ensure monotonicity. 
In our case we are using a robust mixture which goes by applying the well known 
minmod rule to A J+1 / 2 and 2Aj_ 1 / 2 . In that way, the limiter acts only in steep 
regions, where the proportion between neighbouring increments exceeds a factor 
of two. 

We can apply this method to the step-like initial data propagating along 
the x axis. We can see in Figure 3 that the result is much better than before. 
It can be (hardly) observed a small deviation from the TVD condition, which 
is produced by the artificial separation produced by the Strang splitting into 
transport and non-linear source steps. 

This method can be applied, with the general decomposition described in 
section 4, to discontinuities which propagate along any constant direction, and 
not only to the trivial case, aligned with the x axis, that we have considered 
until now. To prove it, we have rotated the metric of Minkowski waves in the 
x-z plane to obtain a diagonal propagation of the profile. The line element in 
this case has the following form: 



ds 2 



x + z 



t) dt 2 + + H{ 



x + z 



t)] (dx 2 +dz 2 )+ dy 2 



V2 



V2 




Figure 4: 3D plot of K xx . The step-like profile has been propagated with 
periodic boundary conditions until one full period (about 80 iterations in this 
case) has elapsed and it has returned to the initial position 

+ h-l + H(?+± - t)] (dx dz + dz dx) (31) 
2 y2 

We show the results in the Figure 4. We can also see in Figure 5 a z=constant 
section of the same results to allow a more detailed comparison with the initial 
data. 
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